Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation.

Undertake the analysis of the primary outcome for ASCOT.

Author

James Totterdell

Published

May 24, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)

theme_set(theme_minimal(base_size = 12))
bayesplot_theme_set(theme_minimal(base_size = 12))
color_scheme_set("red")
options(digits = 4)
Prepare dataset
source("r/derive_full_datasets.r")
all_dat <- read_all_no_daily()

Primary Outcome Definition

The primary outcome is a composite of death, need for new respiratory support, or vasopressor/inotropic support. From the protocol:

Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation. This includes any participant who receives non-invasive mechanical ventilation (either CPAP or BIPAP, apart from the below considerations) any time after enrolment even if not transferred to ICU. It does NOT include the use of humidified high-flow nasal prong oxygen.

Participants on pre-existing home BiPAP or CPAP will not be considered to have met the primary outcome unless they have either:

  • required invasive mechanical ventilation (i.e. intubation), or
  • graduated from CPAP only whilst asleep to BiPAP at any time, or
  • graduated from BiPAP only whilst asleep to BiPAP for >12 hours/day, or
  • died by day 28

There may be cases where a patient has been assessed as requiring intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support, but the patient or family declined treatment and the patient was discharged home. If attempts to obtain 28-day data are unsuccessful or not possible, and the investigator had deemed at the time of discharge that the patient would be highly likely to die within 28 days from randomisation, these participants will be deemed to have met the primary outcome.

Derivation of the Outcome

Derivation of the outcome requires checking of the daily, discharge, and day 28 data extracts. On the daily data, there is a variable DD_PrimaryEndpointReachedToday however this was coded incorrectly in the original database and therefore fails to capture some participants. Additionally, given the composite nature of the outcome it is useful to check all components individually as well as the composite outcome. Therefore, this variable is not used in the derivation, but is included as a cross-check.

Below, each component is summarised in aggregate.

Day 28 mortality

There were 50 participants who had died by day 28, 40 were are reported as a discharge status of death and 10 were reported as a death post-discharge. Day 28 mortality was missing or unknown for 30 participants.

Summarise the death component
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  select(DIS_death, DIS_day, D28_death) %>%
  dplyr::count(DIS_death, DIS_day <= 28, D28_death) %>%
  spread(D28_death, n, fill = 0)
# A tibble: 4 × 5
  DIS_death `DIS_day <= 28`   `0`   `1` `<NA>`
      <int> <lgl>           <dbl> <dbl>  <dbl>
1         0 FALSE              10     0      0
2         0 TRUE             1489    10     30
3         1 FALSE               3     0      0
4         1 TRUE                0    40      0

Vasopressor/inotropic requirement

Vasopressor/inotropic support is reported on both the daily (from day 1 to discharge) and the day 28 (from discharge to day 28) extracts. There were 39 participants who required vasopressor/inotropic support. Of these, 19 were reported between discharge and day 28, 18 were reported prior to discharge, and 2 were reported for both. The component was missing for 35 participants. Either due to missing information at day 28 (33) or due to missing daily information (2).

Summarise the vasopressor/inotropes component
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  select(DAILY_missing, ANY_vasop, ANY_DD_vasop, D28_vasop) %>%
  dplyr::count(DAILY_missing, ANY_DD_vasop, D28_vasop, ANY_vasop) %>%
  spread(ANY_vasop, n, fill = 0)
# A tibble: 7 × 6
  DAILY_missing ANY_DD_vasop D28_vasop   `0`   `1` `<NA>`
          <dbl>        <dbl>     <dbl> <dbl> <dbl>  <dbl>
1             0            0         0  1508     0      0
2             0            0         1     0    19      0
3             0            0        NA     0     0     33
4             0            1         0     0    18      0
5             0            1         1     0     2      0
6             1           NA         0     0     0      1
7             1           NA        NA     0     0      1

New intensive respiratory support

Respiratory support is reported on both the daily (from day 1 to discharge) and the day 28 (at day 28) extracts. There were 70 participants who required new intensive respiratory support. Of these, 60 were prior to discharge, 6 were reported at day 28, and 4 were both. The outcome was unknown for 36 participants due to either missing day 28 information (34) or missing daily information (2).

Summarise the ventilation component
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  select(DAILY_missing, ANY_vent, ANY_DD_vent, D28_vent) %>%
  dplyr::count(DAILY_missing, ANY_DD_vent, D28_vent, ANY_vent) %>%
  spread(ANY_vent, n, fill = 0)
# A tibble: 6 × 6
  DAILY_missing ANY_DD_vent D28_vent   `0`   `1` `<NA>`
          <dbl>       <dbl>    <dbl> <dbl> <dbl>  <dbl>
1             0           0        0  1476     0      0
2             0           0       NA     0     0     34
3             0           1        0     0    60      0
4             0           1        1     0     4      0
5             0           1       NA     0     6      0
6             1          NA        0     0     0      2

DAMA and Likely to Die

There were 3 participants who were DAMA and identified as likely to die. However, the day 28 status was observed for all 3 participants, therefore, no-one met the primary endpoint due to DAMA and unknown day 28 status.

Summarise the DAMA component
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  select(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
  dplyr::count(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
  spread(D28_death, n)
# A tibble: 3 × 5
  DIS_DAMA DIS_DAMAlikelytodie   `0`   `1` `<NA>`
     <int>               <int> <int> <int>  <int>
1        0                   0  1469    46     28
2        1                   0    32     2      2
3        1                   1     1     2     NA

Composite

There were 107 (7% of observed outcomes) participants who met the composite outcome and 36 (2% of all outcomes) participants whose outcome was unknown. Of the 107 who met the outcome, 12 met all three criteria (death, vasopressor, ventilation) and the remaining met some subset of them. Of the 36 with missing outcome, 26 were missing all three components.

Summarise the composite
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(PO)
# A tibble: 3 × 2
     PO     n
  <dbl> <int>
1     0  1439
2     1   107
3    NA    36
Summarise the composite
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  select(PO, D28_death, ANY_vasop, ANY_vent) %>%
  dplyr::count(D28_death, ANY_vasop, ANY_vent, PO) %>%
  spread(PO, n, fill = 0)
# A tibble: 15 × 6
   D28_death ANY_vasop ANY_vent   `0`   `1` `<NA>`
       <int>     <dbl>    <dbl> <dbl> <dbl>  <dbl>
 1         0         0        0  1439     0      0
 2         0         0        1     0    31      0
 3         0         1        0     0    20      0
 4         0         1        1     0     5      0
 5         0        NA        0     0     0      5
 6         0        NA       NA     0     0      2
 7         1         0        0     0    10      0
 8         1         0        1     0    21      0
 9         1         0       NA     0     4      0
10         1         1        0     0     2      0
11         1         1        1     0    12      0
12         1        NA       NA     0     1      0
13        NA         0       NA     0     0      3
14        NA        NA        1     0     1      0
15        NA        NA       NA     0     0     26

Descriptive Analyses

Age

Relationship age to outcome
agedat <- all_dat %>% 
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(PO, AgeAtEntry) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion with outcome", x = "Age at entry")
ggplotly(p)

Relationship age 60+ to outcome
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(`Age 60+` = agegte60, PO) %>%
  group_by(`Age 60+`) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped")
Table 1: Relationship between age >= 60 and the primary outcome.
Age 60+ 0 1 <NA> p_1 p_miss
0 1038 59 25 0.05 0.02
1 401 48 11 0.11 0.02

Country

Relationship country to outcome
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(Country = PT_CountryName, PO) %>%
  group_by(Country) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped")
Table 2: Primary outcome by country.
Country 0 1 <NA> p_1 p_miss
Australia 128 8 14 0.06 0.09
India 1190 83 15 0.07 0.01
Nepal 104 11 3 0.10 0.03
New Zealand 17 5 4 0.23 0.15

Site

Relationship site to outcome
tabdat <- all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(
    Country = PT_CountryName, 
    Site = PT_LocationName, 
    PO) %>%
  group_by(Country, Site) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()
n_country <- tabdat %>% 
  dplyr::count(Country) %>%
  pull(n)
cn <- cumsum(n_country)
tabdat %>%
  select(-Country) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped") %>%
  pack_rows("Australia", 1, cn[1]) %>%
  pack_rows("India", cn[1] + 1, cn[2]) %>%
  pack_rows("Nepal", cn[2] + 1, cn[3]) %>%
  pack_rows("New Zealand", cn[3] + 1, cn[4])
Table 3: Primary outcome by site within country.
Site 0 1 <NA> p_1 p_miss
Australia
Alfred Hospital 6 1 7 0.14 0.50
Blacktown Hospital 4 0 0 0.00 0.00
Box Hill Hospital 4 2 2 0.33 0.25
Campbelltown Hospital 11 1 0 0.08 0.00
John Hunter Hospital 2 0 0 0.00 0.00
Liverpool Hospital 6 2 0 0.25 0.00
Monash Health 8 0 0 0.00 0.00
Prince Charles Hospital 3 0 0 0.00 0.00
Royal Melbourne Hospital 7 1 1 0.12 0.11
Royal North Shore Hospital 6 0 0 0.00 0.00
Royal Perth Hospital 5 0 1 0.00 0.17
St George Hospital 1 0 0 0.00 0.00
St Vincent's Hospital, Sydney 10 0 1 0.00 0.09
Wagga Wagga Base Hospital 10 0 0 0.00 0.00
Western Health 1 0 0 0.00 0.00
Westmead Hospital 42 1 2 0.02 0.04
Wollongong Hospital 2 0 0 0.00 0.00
India
Aditya Multispeciality Hospital 63 0 0 0.00 0.00
Apex Hospitals, Malviya Nagar 125 1 0 0.01 0.00
Believers Church Medical College Hospital 18 0 0 0.00 0.00
Christian Medical College Vellore 84 16 10 0.16 0.09
CMC Ludhiana 64 18 0 0.22 0.00
Core Hospital 86 21 0 0.20 0.00
Jivanrekha Multispeciality Hospital 329 13 3 0.04 0.01
Maharaja Agrasen Superspeciality Hospital 138 0 0 0.00 0.00
Samishta Hospital & Research Institute 156 2 0 0.01 0.00
Sterling Multispecialty Hospital Pune 85 12 2 0.12 0.02
Symbiosis University Hospital 42 0 0 0.00 0.00
Nepal
Bir Hospital 71 8 3 0.10 0.04
Tribhuvan University Teaching Hospital 33 3 0 0.08 0.00
New Zealand
Christchurch Hospital 1 0 0 0.00 0.00
Dunedin Hospital 1 0 1 0.00 0.50
Middlemore Hospital 13 5 0 0.28 0.00
North Shore Hospital 2 0 3 0.00 0.60

Calendar Time

Relationship calendar date to outcome
caldat <- all_dat %>% 
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(PO, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`))
p <- ggplot(caldat, aes(mth, p)) +
  facet_wrap( ~ yr) +
    geom_point() +
    labs(
      y = "Proportion with outcome", 
      x = "Calendar date (year and month)") +
  scale_x_continuous(breaks = 1:12)
ggplotly(p)

Anti-coagulation Domain

Relationship anti-coagulation to outcome
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(`Anti-coagulation` = CAssignment, PO) %>%
  group_by(`Anti-coagulation`) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped")
Table 4: Primary outcome by anti-coagulation intervention.
Anti-coagulation 0 1 <NA> p_1 p_miss
C0 16 4 6 0.20 0.23
C1 551 45 14 0.08 0.02
C2 571 30 12 0.05 0.02
C3 258 21 4 0.08 0.01
C4 43 7 0 0.14 0.00
Relationship anti-coagulation to outcome by country
all_dat %>%
  filter(ENR_rec == 1, WTH_FU == 0) %>%
  dplyr::count(
    Country = PT_CountryName, 
    `Anti-coagulation` = CAssignment, 
    PO) %>%
  complete(Country, `Anti-coagulation`, PO, fill = list(n = 0)) %>%
  group_by(Country, `Anti-coagulation`) %>%
  spread(PO, n, fill = 0) %>%
  ungroup() %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  select(-Country) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped") %>%
  pack_rows("Australia", 1, 5) %>%
  pack_rows("India", 6, 10) %>%
  pack_rows("Nepal", 11, 15) %>%
  pack_rows("New Zealand", 16, 20)
Table 5: Primary outcome by anti-coagulation intervention by country.
Anti-coagulation 0 1 <NA> p_1 p_miss
Australia
C0 16 2 6 0.11 0.25
C1 45 2 2 0.04 0.04
C2 50 3 6 0.06 0.10
C3 7 0 0 0.00 0.00
C4 10 1 0 0.09 0.00
India
C0 0 0 0 NaN NaN
C1 449 37 7 0.08 0.01
C2 486 25 5 0.05 0.01
C3 251 21 3 0.08 0.01
C4 4 0 0 0.00 0.00
Nepal
C0 0 0 0 NaN NaN
C1 49 4 3 0.08 0.05
C2 30 1 0 0.03 0.00
C3 0 0 0 NaN NaN
C4 25 6 0 0.19 0.00
New Zealand
C0 0 2 0 1.00 0.00
C1 8 2 2 0.20 0.17
C2 5 1 1 0.17 0.14
C3 0 0 1 NaN 1.00
C4 4 0 0 0.00 0.00

Analyses

In what follows, the model is built up from the simplest (treatment-only) to the most complex (primary model).

ACS-ITT

Treatment only

Data and sampling treatment only model
seed <- 17250
logistic_mod <- cmdstan_model("stan/logistic_bernoulli.stan")

make_model_data_trt_only <- function(dat) {
  X <- model.matrix(
     ~ CAssignment,
    data = dat[!is.na(dat$PO), ],
    contrast = list(CAssignment = contr.orthonorm)
  )
  y <- as.integer(dat$PO)
  y <- y[!is.na(y)]
  N <- dim(X)[1]
  K <- dim(X)[2]  
  list(N = N, K = K, X = X, y = y, beta_sd = c(2.5, 1, 1, 1))
}
mod_dat <- make_model_data_trt_only(
  all_dat %>%
    filter(ENR_rec == 1, ACS_ITT == 1, WTH_FU == 0)
)

snk <- capture.output(
  trt_only_fit <- logistic_mod$sample(
    data = mod_dat,
    seed = seed,
    refresh = 0,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8
))

Diagnostics

     num_divergent num_max_treedepth ebfmi
[1,]             0                 0 1.073
[2,]             0                 0 1.088
[3,]             0                 0 1.029
[4,]             0                 0 1.125
[5,]             0                 0 1.000
[6,]             0                 0 1.054
[7,]             0                 0 1.103
[8,]             0                 0 1.096
Summaries
trt_only_fit$summary("beta")
# A tibble: 4 × 10
  variable    mean  median    sd   mad       q5     q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -2.48   -2.48   0.132 0.131 -2.70    -2.27    1.00   12751.   13864.
2 beta[2]   0.574   0.582  0.344 0.343 -0.00605  1.12    1.00   13362.   13385.
3 beta[3]  -0.0352 -0.0360 0.196 0.198 -0.357    0.287   1.00   13470.   14218.
4 beta[4]  -0.413  -0.413  0.209 0.209 -0.756   -0.0712  1.00   14124.   13591.
Trace plots
mcmc_trace(trt_only_fit$draws("beta_raw", format = "matrix"))

Posterior draws
C <- cbind(1, contr.orthonorm(4))
modrv <- as_draws_rvars(trt_only_fit$draws("beta"))
modrv$eta <- as.vector(C %**% modrv$beta[1:4])
modrv$beta_uncon <- as.vector(C[, -1] %**% modrv$beta[2:4])
modrv$prob <- 1 / (1 + exp(-modrv$eta))
modrv$beta_trt <- as.vector(
  transform_coding(C, cbind(1, contr.treatment(4))) %**% 
  modrv$beta[1:4])
modrv$OR <- setNames(exp(modrv$beta_trt[-1]), c("C2", "C3", "C4"))
Posterior density
mcmc_hist(modrv["OR"], binwidth = 0.02) +
  labs(x = "Odds ratio (vs C1)") +
  geom_vline(xintercept = 1, linetype = 2) +
  scale_x_log10()

Treatment plus covariates (no site or epoch)

Under this model a component is added for country (region), \(\mathbb{I}[\text{age}\geq 60]\) (sample-mean centred), and intervention specific ineligibility (C3 only)

New terms:

\[ \begin{aligned} \rho_1 &= 0 \\ \rho_r &\overset{\text{iid}}{\sim} \text{Normal}(0, 1), r=2,...,R,\\ \beta_{\text{age}\geq 60} &\sim \text{Normal}(0, 2.5^2) \\ \beta_{\text{inelgC3}} &\sim \text{Normal}(0, 10) \end{aligned} \]

Data and sampling no site or epoch model
seed <- 71236
logistic_mod <- cmdstan_model("stan/logistic_bernoulli.stan")

make_model_data_noepochsite <- function(dat) {
  X <- model.matrix(
     ~ CAssignment + agegte60_c + country + inelgc3,
    data = dat[!is.na(dat$PO), ],
    contrast = list(CAssignment = contr.orthonorm)
  )
  y <- as.integer(dat$PO[!is.na(dat$PO)])
  N <- dim(X)[1]
  K <- dim(X)[2]  
  list(N = N, K = K, X = X, y = y, 
       beta_sd = c(2.5, 1, 1, 1, 2.5, 1, 1, 1, 10))
}
mod_dat <- make_model_data_noepochsite(
  all_dat %>%
    filter(ENR_rec == 1, ACS_ITT == 1, WTH_FU == 0) %>%
    mutate(
      age_c = AgeAtEntry - mean(AgeAtEntry),
      agegte60_c = agegte60 - mean(agegte60),
      country = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
      inelgc3 = if_else(EL_inelg_c3 == 1 | is.na(EL_inelg_c3), 1, 0)
    )
)

snk <- capture.output(
  no_epoch_site_fit <- logistic_mod$sample(
    data = mod_dat,
    seed = seed,
    refresh = 0,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8
))
Summaries
no_epoch_site_fit$summary("beta")
# A tibble: 9 × 10
  variable    mean  median    sd   mad     q5     q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -2.49   -2.48   0.181 0.180 -2.79  -2.19    1.00   12045.   13850.
2 beta[2]   0.547   0.550  0.403 0.402 -0.122  1.20    1.00   13742.   13325.
3 beta[3]   0.0341  0.0338 0.202 0.202 -0.294  0.370   1.00   16197.   15653.
4 beta[4]  -0.323  -0.324  0.212 0.211 -0.669  0.0217  1.00   21503.   16268.
5 beta[5]   0.760   0.759  0.209 0.209  0.417  1.10    1.00   22417.   15132.
6 beta[6]  -0.282  -0.269  0.407 0.408 -0.974  0.358   1.00   21069.   14367.
7 beta[7]   0.248   0.253  0.391 0.388 -0.402  0.883   1.00   16499.   15453.
8 beta[8]   0.573   0.601  0.597 0.587 -0.447  1.51    1.00   21254.   13934.
9 beta[9]  -0.304  -0.297  0.260 0.258 -0.742  0.116   1.00   16008.   15128.

Diagnostics

     num_divergent num_max_treedepth  ebfmi
[1,]             0                 0 1.0572
[2,]             0                 0 1.0200
[3,]             0                 0 1.0122
[4,]             0                 0 0.9957
[5,]             0                 0 1.1094
[6,]             0                 0 1.0026
[7,]             0                 0 1.0076
[8,]             0                 0 1.0577
Trace plots
mcmc_trace(no_epoch_site_fit$draws("beta_raw", format = "matrix"))

Posterior draws
C <- cbind(1, contr.orthonorm(4))
modrv <- as_draws_rvars(no_epoch_site_fit$draws("beta"))
modrv$eta <- as.vector(C %**% modrv$beta[1:4])
modrv$beta_uncon <- as.vector(C[, -1] %**% modrv$beta[2:4])
modrv$prob <- 1 / (1 + exp(-modrv$eta))
modrv$beta_trt <- as.vector(
  transform_coding(C, cbind(1, contr.treatment(4))) %**% 
  modrv$beta[1:4])
modrv$OR <- setNames(exp(modrv$beta_trt[-1]), c("C2", "C3", "C4"))
Posterior density
mcmc_hist(modrv["OR"], binwidth = 0.02) +
  labs(x = "Odds ratio (vs C1)") +
  geom_vline(xintercept = 1, linetype = 2) +
  scale_x_log10()

Posterior summaries
tibble(
    Parameter = c("OR[C2]", "OR[C3]", "OR[C4]"),
    `Pr(OR < 1)` = Pr(modrv$OR < 1),
    Median = median(modrv$OR),
    Mean = E(modrv$OR),
    `95% CrI` = apply(
      quantile(modrv$OR, c(0.025, 0.975)), 2, 
      function(x) sprintf("(%.2f, %.2f)", x[1], x[2]))
  ) %>%
  kable(digits = 2, align = "lrrrr") %>%
  kable_styling(bootstrap_options = "striped")
Table 6: Posterior summaries
Parameter Pr(OR < 1) Median Mean 95% CrI
OR[C2] 0.96 0.65 0.67 (0.40, 1.04)
OR[C3] 0.68 0.88 0.91 (0.50, 1.53)
OR[C4] 0.17 1.58 1.75 (0.60, 3.92)

Treatment plus covariates plus site (no epoch)

Under this model a component is added for site nested within country.

New terms

\[ \begin{aligned} \xi_{rs} &\overset{\text{iid}}{\sim}\text{Normal}(0,\sigma^2_{r}), s=1,...,S_r, r=1,...,R \\ \sigma_r &\overset{\text{iid}}{\sim} \text{Half-}t(3, 1). \end{aligned} \]

Data and sampling no site or epoch model
seed <- 71236
noepoch_mod <- cmdstan_model("stan/primary_logistic_model_site.stan")
make_model_data_noepoch <- function(dat) {
  
  site_counts <- dat %>% 
    dplyr::count(
      Country = factor(Country, levels = c("IN", "AU", "NP", "NZ")), 
      Location)
  merge_aus <- site_counts %>% 
    filter(Country == "AU", n < 5) %>%
    pull(Location)
  merge_nz <- site_counts %>% 
    filter(Country == "NZ", n < 5) %>%
    pull(Location)
  dat <- dat %>%
    mutate(
      ctry = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
      # group sites with < 5 counts
      site = fct_collapse(
        factor(Location, levels = site_counts$Location),
        `AU other` = merge_aus,
        `NZ other` = merge_nz
        )
    )
  region_site <- dat %>%
      dplyr::count(ctry, site) %>%
      mutate(
        ctry_num = as.numeric(ctry), 
        site_num = as.numeric(fct_inorder(site))
      )
  dat <- dat %>% 
    left_join(region_site %>% select(-n), by = c("ctry", "site"))
  region <- dat$ctry_num
  M_region <- max(region)
  site <- dat$site_num
  M_site <- max(site)
  region_by_site <- region_site$ctry_num
  
  X <- model.matrix(
     ~ CAssignment + agegte60_c + inelgc3,
    data = dat,
    contrast = list(CAssignment = contr.orthonorm)
  )
  y <- as.integer(dat$PO)
  N <- dim(X)[1]
  K <- dim(X)[2]  
  list(N = N, K = K, X = X, y = y,
       M_region = M_region, region = region,
       M_site = M_site, site = site,
       region_by_site = region_by_site,
       beta_sd = c(2.5, 1, 1, 1, 2.5, 10))
}

sub_dat <- all_dat %>%
    filter(
      ENR_rec == 1, 
      ACS_ITT == 1, 
      WTH_FU == 0,
      !is.na(PO)) %>%
    mutate(
      age_c = AgeAtEntry - mean(AgeAtEntry),
      agegte60_c = agegte60 - mean(agegte60),
      country = factor(Country, levels = c("IN", "AU", "NP", "NZ")),
      inelgc3 = if_else(EL_inelg_c3 == 1 | is.na(EL_inelg_c3), 1, 0)
    )

mod_dat <- make_model_data_noepoch(sub_dat)

snk <- capture.output(
  no_epoch_fit <- noepoch_mod$sample(
    data = mod_dat,
    seed = seed,
    refresh = 0,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8,
    adapt_delta = 0.9
))

Diagnostics

     num_divergent num_max_treedepth  ebfmi
[1,]             0                 0 0.8572
[2,]             0                 0 0.7903
[3,]             0                 0 0.8450
[4,]             0                 0 0.7790
[5,]             0                 0 0.8532
[6,]             0                 0 0.7963
[7,]             0                 0 0.7857
[8,]             0                 0 0.8210
Summaries
no_epoch_fit$summary("beta")
# A tibble: 6 × 10
  variable   mean median    sd   mad      q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -3.37  -3.36  0.521 0.515 -4.23   -2.54   1.00    5013.    7718.
2 beta[2]   0.680  0.684 0.428 0.430 -0.0357  1.38   1.00   19585.   15465.
3 beta[3]   0.106  0.105 0.214 0.216 -0.245   0.457  1.00   18640.   15799.
4 beta[4]  -0.258 -0.256 0.228 0.228 -0.631   0.116  1.00   22127.   14602.
5 beta[5]   0.682  0.684 0.229 0.231  0.306   1.06   1.00   25053.   14777.
6 beta[6]   0.478  0.483 0.309 0.308 -0.0355  0.978  1.00   18180.   15850.
Trace plots
mcmc_trace(no_epoch_fit$draws(
  c("beta_raw", "beta_region_raw"), 
  format = "matrix"))

Trace plots site terms
mcmc_trace(no_epoch_fit$draws(
  c("epsilon_site", "tau_site"), 
  format = "matrix"))

Posterior draws
C <- cbind(1, contr.orthonorm(4))
modrv <- as_draws_rvars(
  no_epoch_fit$draws(c("beta", "beta_region", "gamma_site", "tau_site")))
modrv$eta <- as.vector(C %**% modrv$beta[1:4])
modrv$beta_uncon <- as.vector(C[, -1] %**% modrv$beta[2:4])
modrv$prob <- 1 / (1 + exp(-modrv$eta))
modrv$beta_trt <- as.vector(
  transform_coding(C, cbind(1, contr.treatment(4))) %**% 
  modrv$beta[1:4])
modrv$OR <- setNames(exp(modrv$beta_trt[-1]), c("C2", "C3", "C4"))
Posterior density
mcmc_hist(modrv["OR"], binwidth = 0.02, freq = F) +
  labs(x = "Odds ratio (vs C1)") +
  geom_vline(xintercept = 1, linetype = 2) +
  scale_x_log10()

Posterior summaries
tibble(
    Parameter = c("OR[C2]", "OR[C3]", "OR[C4]"),
    `Pr(OR < 1)` = Pr(modrv$OR < 1),
    Median = median(modrv$OR),
    Mean = E(modrv$OR),
    `95% CrI` = apply(
      quantile(modrv$OR, c(0.025, 0.975)), 2, 
      function(x) sprintf("(%.2f, %.2f)", x[1], x[2]))
  ) %>%
  kable(digits = 2, align = "lrrrr") %>%
  kable_styling(bootstrap_options = "striped")
Table 7: Posterior summaries
Parameter Pr(OR < 1) Median Mean 95% CrI
OR[C2] 0.98 0.61 0.62 (0.37, 0.99)
OR[C3] 0.85 0.72 0.76 (0.40, 1.30)
OR[C4] 0.18 1.60 1.81 (0.57, 4.21)

Primary Model

FAS-ITT